arXiv:1509.02015v2 [math.NA] 18 Mar 2016 


METHODS FOR VERIFIED STABILIZING SOLUTIONS TO 
CONTINUOUS-TIME ALGEBRAIC RICCATI EQUATIONS 

TAYYEBE HAQIRI AND FEDERICO POLONI 


Abstract. We describe a procedure based on the Krawczyk method to 
compute a verified enclosure for the stabilizing solution of a continuous¬ 
time algebraic Riccati equation A*X + XA + Q — XGX building on the 
work of [B. Hashemi, SCAN 2012] and adding several modifications to 
the Krawczyk procedure. We show that after these improvements the 
Krawczyk method reaches results comparable with the current state-of-the- 
art algorithm [Miyajima, Jpn. J. Ind. Appl. Math 2015], and surpasses 
it in some examples. Moreover, we introduce a new direct method for 
verification which has a cubic complexity in term of the dimension of X, 
employing a fixed-point formulation of the equation inspired by the ADI 
procedure. The resulting methods are tested on a number of standard 
benchmark examples. 


[2010]65M32, 35Kxx, 65T60. Algebraic Riccati equation, stabilizing solu¬ 
tion, interval arithmetic, verified computation, Krawczyk’s method. 


1. Introduction 

Consider the continuous-time algebraic Riccati equation CARE 

A*X + XA + Q = XGX, (1.1) 

where A, G and Q € are given, G and Q are Hermitian, and X € is 
unknown. Here, the notation A* denotes the conjugate transpose of a complex 
matrix A while shows the transpose of A. CAREs have a variety of applica¬ 
tions in the field of control theory and filter design, such as the linear-quadratic 
optimal control problem and Hamiltonian systems of differential equations [20] . 
A solution Xg of (1.1) is called stabilizing if the closed loop matrix A — GXg 
is Hurwitz stable, i.e., if all its eigenvalues have strictly negative real part. If 
a stabilizing solution Xg exists, it is unique [5, Theorem 2.17], and it is Her¬ 
mitian, i.e., Xg = (Xg)*. The unique stabilizing solution is the one of interest 
in almost all applications [20, 5], hence in this paper we focus on its computa¬ 
tion. The techniques presented here can be adapted with minor sign changes 
to anti-stabilizing solutions, i.e., solutions Xag for which all the eigenvalues of 
A — GXas have positive real part. The algorithms in [14] and [25], in contrast. 
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do not restrict to verifying stabilizing solutions only; however, solutions which 
are neither stabilizing nor anti-stabilizing have very few applicative uses. 

The solutions of (1-1) can be put in one-to-one correspondence with certain 
invariant subspaces of the Hamiltonian 


H : = 


A 

-Q 


-G 

-A* 


€ C 


2nx2n 


Indeed, X is a solution of (1.1) if and only if 


'In 


'In 

X 


X 


{A-GX), 


( 1 . 2 ) 


in which is the identity matrix in In particular, the columns of the 

matrix [ ^ ] span an invariant subspace for the matrix H, and the eigenvalues 
oi A — GX are a subset of the eigenvalues of H [5] . We refer the reader to the 
books by Lancaster and Rodman [20] and by Bini, lannazzo and Meini [5] for 
details concerning main theoretical properties and numerical solutions together 
with the description of the main tools for the design and analysis of solution 
algorithms. 

The work presented in this paper addresses the problem of computing ver¬ 
ified stabilizing solutions of CAREs (1.1), that is, determining an interval 
matrix which is guaranteed to contain the stabilizing solution of the CARE. 
The problem of computing verified solutions to the matrix Riccati equations 
(AREs) has been addressed before in the literature: the algorithms in [21] 
and [22], based on the interval Newton method, are pioneering works in this 
context but their computational complexity is 0{n^). In [21], the authors ap¬ 
ply Brouwer’s fixed point theorem to calculate verified solutions of the ARE 


A^X + XA + Q = XBR-^B'^X, (1.3) 

with real symmetric matrices Q and R, Q positive semi dehnite and R posi¬ 
tive definite. They find an interval matrix including a positive definite solution 
of (1.3). The paper [31] decreases this cost to 0(n^) by using the Krawczyk 
method, which is a variant of the Newton method that does not require the 
inversion of an interval matrix. A major improvement is the algorithm in [14], 
which is applicable when the closed-loop matrix A — GX is diagonalizable 
where X denotes a numerical computed solution of (1.1), and requires only 
0{n^) operations. The recent paper [25] describes a more efficient algorithm 
based again on the diagonalization of A — GX, X a Hermitian numerical so¬ 
lution of (1.1). The resulting method has cubic complexity as well. An im¬ 
portant feature of this algorithm is that does not require iteration to hnd a 
suitable candidate interval solution, unlike the previous methods based on the 
Krawczyk method and hxed-point theorems. Hence it is typically faster than 
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the alternatives. The same paper [25] also includes a method to verify the 
uniqueness and the stabilizing property of the computed solution. 

We propose here a variant of the Krawczyk method suggested in [14] , intro¬ 
ducing several modifications. Namely: 

• We use the technique introduced in [9] , which consists in applying the 
Krawczyk method not to the original equation, but to one obtained 
after a change of basis, in order to reduce the number of verified oper¬ 
ations required, with the aim to reduce the wrapping effects. 

• We exploit the invariant subspace formulation (1.2) to make another 
change of basis in the matrix H, following a technique introduced 
in [23] for the non-verified solution of Riccati equations. This tech¬ 
nique employs suitable permutations of H to reduce (1.1) to a differ¬ 
ent CARE whose stabilizing solution Yg has bounded norm. Up to our 
knowledge, this is the first attempt to use this technique in the context 
of interval arithmetic and verified computation. 

• When applying the Krawczyk method, an enclosure for the so-called 
slope matrix is needed; the standard choice to compute it is using the 
interval evaluation of the Jacobian of the function at hand. Instead, we 
use a different algebraic expression which results in a smaller interval. 

An algorithm to verify the stabilizing property (and thus the uniqueness) of 
the solution enclosure in the computed interval matrix is also presented. 

In addition, we present a different algorithm, based on a reformulation 
of (1.1) as a fixed-point equation, which requires 0{n^) operations per step 
and does not require the diagonalizability of the closed-loop matrix A — GX 
in which X is the computed approximate stabilizing solution of CARE (1.1). 
This algorithm is generally less reliable than the Krawczyk-based ones, but 
it has the advantage of not breaking down in cases in which the closed-loop 
matrix is defective or almost defective. 

We conclude the paper by evaluating the proposed algorithms on a large set 
of standard benchmark problems [4, 7] for Riccati equations, comparing them 
with the algorithms in [14] and [25]. Using all the improvements described here, 
the gap between the Krawczyk method and the current best method in [25] is 
essentially eliminated. The four methods each handle satisfactorily a slightly 
different set of problems, and none of them is beaten by the alternatives in all 
possible experiments. 

The paper is organized as follows. In the next section we introduce some 
notation and standard results in linear algebra and interval analysis which 
are at the basis of our methods. In section 3 we discuss various algorithms 
based on the Krawczyk method to compute a thin interval matrix enclosing 
a solution of (1.1) while in section 4, a fixed point approach is presented. In 
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section 5 and 6 we perform some numerical tests and draw the conclusions 
and outlook, respectively. 

2. Preliminaries and Notation 

We try to follow the standard notation of interval analysis defined in [18]. 
Subsequently, we use boldface lower and upper case letters for interval scalars 
or vectors and matrices, respectively, whereas lower case stands for scalar quan¬ 
tities and point vectors and upper case represents matrices. By 
and we denote the sets of all complex (real) m x n matrices 

and the set of all complex (real) m x n interval matrices, respectively. By a 
complex (real) interval matrix, we mean a matrix whose entries are circular 
complex (compact real) intervals. 

The Kronecker product A® B of an m x n matrix A = {Aij) and a p x q 
matrix B is an mp x nq matrix defined as the block matrix whose blocks are 
A® B ■.= [AijB], For a point matrix A G the vector vec(74) € C™” 

denotes column-wise vectorization whereby the successive columns of A are 
stacked one below the other, beginning with the first column and ending with 
the last. Moreover, A denotes the complex conjugate of A and if A is an 
invertible matrix, then A~'^ := and A~* := (^*)“^. The element-wise 

division of a matrix A = {Aij) € c^xn |^y matrix B = {Bij) G 
also known as the Hadamard division, denoted by A.jB, results in an m x n 
matrix C = {Cij) whose (i,j)-th element is given by Cij = AijjBij provided 
that Bij 7 ^ 0, for each 1 < i < m and 1 < j < n. For a given vector 
d = {di,d 2 , ... ,dn)'^ G C"’, diag((i) G is the diagonal matrix whose 

entry is di. Conversely, given a diagonal matrix D, diag(T>) is the vector whose 
elements are the diagonal entries of D. Most of these notions and operations 
are analogously defined for interval quantities. 

Complex intervals can be defined either as rectangles or as discs. We use 
here the definition as discs: a circular complex interval x, or circular disc 
or simply a complex interval, is a closed circular disc of radius rad (x) G M 
with rad (x) > 0 and center mid (x) G C. Indeed, it is defined as x := {z G 
C: |z—mid(x)| <rad(x)} = (mid (x), rad (x)). The operations on the circular 
complex intervals, ICdiso ^-^e introduced as generalizations of operations on 
complex numbers. Then, the standard arithmetic for circular complex interval 
arguments x = (mid (x), rad (x)) and y = (mid (y), rad (y)) is defined as [2]: 

X ± y := (mid (x) ± mid (y), rad (x) -|- rad (y)) , 

x-y := (mid (x) mid (y), | mid (x)| rad (y) -|- | mid (y)| rad (x) -|- rad (x) rad (y)) , 
1/x := ^mid (x)/| mid (x)p — (rad (x))^,rad(x)/| mid (x)|^ — (rad (x))^^ , 0 ^ x, 

x/y:=x-l/y, 0 ^ y. 
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Here \z\ denotes the absolute value of a complex number z. Additionally, op¬ 
erations between a complex interval and a complex number z G C can be 
performed by identifying z with (z, 0) € ICdisc which can be then extended to 
other point and interval quantities. Besides, z is the complex conjugate of z. 

We emphasize that the definitions of addition, subtraction and inversion sub¬ 
stantially coincide with their set theoretic definitions, i.e., 

X ± y = {x ± y: X G X, y G y}, 

l/x = {l/x:xGx}. 

Unfortunately, the set {xy: x G x, y G y}, in general, is not a disc but one of 
the basic properties of interval arithmetic, which makes its use well-founded, 
is that respects inclusion: for all the four basic arithmetic operations o G 
{-|-, —, •, /} one has 

xoy D {xoy : x G x, y G y}, 

in which x and y are two real or circular complex intervals. (In the case of 
division, we need to assume that 0 ^ y for the operation to be well-defined.) 
Hence we have the following inclusion property: if f is the interval evaluation 
of / (an ordinary real or complex function of N variables), then 

/(xi,X 2 ,... ,XAr) := {/(xi,X 2 ,... ,XAr) : xi G xi,X 2 G X 2 ,... ,XAr G xat} 

C f(xi,X2,... ,XAr). 

The interval evaluation [26, e.g.] of a function /(xi, X 2 ,..., xat), f (xi, X 2 ,..., x^r), 
defined by a formula is obtained by replacing (1) the variables xi, X 2 , • • •, xn 
with interval variables xi, X 2 ,..., xjv and (2) each arithmetic operation with 
the corresponding interval operation. We will utilize the same approach to de¬ 
fine the interval evaluation of a matrix function as well. Note that, in principle, 
different equivalent formulas could give different interval evaluations; indeed, 
the process of turning the customary arithmetic into interval arithmetic is 
not free of pitfalls; issues such as interval dependency and the wrapping phe¬ 
nomenon have to be considered carefully. We refer the reader to the review 
article [28] for a thorough introduction. 

The intersection of two complex intervals x and y is not always a complex 
interval so we may define it to be any complex interval z such that x n y C z, 
if they are not disjoint. If they are disjoint, the intersection is the empty set. 
The function intersect .m in INTLAB can be used to obtain a tight intersect. 
The interval hull of two intervals x and y is denoted by □(x, y) which is the 
smallest interval containing x and y. The magnitude of x G ICdisc is defined 
as mag(x) := max{|x| : x G x}. 

We denote by A = (mid(A),rad(A)) G IC^JIg^"" the m x n interval ma¬ 
trix A whose {i,j) element is the complex interval (mid(Ajj), rad(Ajj)), with 
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rad(Ajj) >0; 1 < i < m,l < j < n. For interval vectors and matrices, 
mid, rad, mag, intersect and □ will be applied component-wise. 

The Frobenius norm of a complex matrix A = (Aij) is defined as ||A||j;. := 

This definition can be extended to complex interval matrices, 
providing an interval-valued function ||A||p defined as the smallest interval 
containing {x : x = ||A||^, A G A}. 

The definition of inverse of an interval matrix may be problematic in general, 
but if D = diag(d) is a diagonal interval matrix, with d = (di, d 2 ,..., d^v)^ 
and 0 0 dj for each i = 1,2,... ,N, then we may define := diag((d]" , d^ , • • •, d^ )^). 

The following lemmas contain simple arithmetical properties of the Kro- 
necker products which we will use in the following. Most of them appear also 
e.g. in [9] or [17]. 

Lemma 2.1. Assume that A = {Aij), B = {Bij), C = {Cij) and D = {Dij) 
he complex matrices with compatible sizes. Then, 

(1) {A^ B){C ^ D) = AC ^ BD, 

(2) A0 {B + C) = (A® B) + {A® C), 

(3) {A® B)* = A* ® B*, 

(4) (A ® B)~^ = A~^ ® if A and B are invertible, 

(5) vec{ABC) = ® A) vec (B), 

(6) (diag(vec (A)))“^ vec (B) = vec (B./A), if Aij / 0 for each i,j. 

Lemma 2.2. Let A = (Aij), B = (Bjj) and C = (Cij) be complex interval 
matrices of compatible sizes. Then, 

( . fvecfA(BC)V 

(1) \{C^ ®A)Yec{B): Ae A, B gB, C eC\ F { ) ( 

^ ^ lvecnAB)Cj, 

(2) (^diag( vec vec (B) = vec (B./A), if 0 ^ Aij for all i,j. 

The CPU rounding mode called round to nearest is typically used when we 
compute using a computer. Two other rounding modes are round downward 
(round towards the largest floating point number smaller than the true result), 
and round upward (round towards the smallest floating point number larger 
than the true result) [1]. In general, it is impossible to compute the exact 
solution in a numerical computation using floating point numbers. Then, a 
numerical computation with guaranteed accuracy switches the CPU rounding 
mode, computes a lower and an upper bound to the true solution, and creates 
an interval which is guaranteed to contain it [13]. One example of software 
which provides a fast implementation of such a reliable interval arithmetic 
is the MATLAB toolbox INTLAB [29]; older versions of INTLAB are freely 
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available for noncommercial use. The default arithmetic for both real and 
complex intervals in INTLAB is the midpoint-radius arithmetic [28]. 

3. Modified Krawczyk’s methods 

Enclosure methods using interval arithmetic work this way: let g : —>■ 

C'^ be some function of which we wish to find a zero. First find a function 
h : —>■ whose fixed points are known to be the zeros of g. Assume 

that h is continuous and that we know an enclosure function h for h, i.e., a 
function based on ‘correct’ interval arithmetic which gives an interval vector 
h(x) containing the range of h over a given interval vector x. Then if h(x) C 
X we know that h(x) C x and so h has a hxed point in x by Brouwer’s 
theorem [8]. 

In this paper, often the functions h are variants of the Krawczyk operator. 
To define this operator, we first need the concept of a slope. 

Definition 3.1. [26, e.g.j Suppose / : V’ ^ —>• and . Then, 

a slope S{f;x,y) is a mapping from the Cartesian product ip x ip to 
such that 

fiy) - fix) = Sif] X, y){y - x). 

We are now ready to state the result which is at the basis of all the modified 
Krawczyk-type algorithms used in the rest of our paper. 

Theorem 3.2. [11] Assume that f : ip G is continuous. Let x ^ pj 

and z G be such that x + z G ip. Moreover, assume that S G is 

a set of matrices such that S{f\x,x') G S for every x' ^ x + z =: x. Finally, 
let R G . Denote by JCf(x,R,z,S) the set 

K,fix,R,z,S) := {—Rf{x) + {In — RS)z : 5 G 5, z G z}. 

If 

ICf{x,R,z,S)Gint{z), (3.1) 

then the function f has a zero x* in x + lCf{x, R,z,S) G x, in which int(z) is 
the topological interior of z. 

Moreover, if S{f;y,y') G S for each y,y' G x, then x* is the only zero of f 
contained in x. 

In computation, one defines the Krawczyk operator [19] 

kj(x, R, z, S) := —Rf{x) + {In — RS)z, (3-2) 

where S is an interval matrix containing all slopes S{f‘,y,y') for y,y' G x. 
A common choice for S is obtained from f^(x), an interval evaluation of the 
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Jacobian f' on the interval x. By the inclusion property of interval arithmetic, 

kj(x, z, S) C int (z) (3.3) 

implies (3.1). So, if (3.3) is satisfied then / has a zero in + kj(x, i?, z, S). In 
practice, one attempts to make the terms —Rf{x) and In — RS as small as 
possible, to obtain the crucial relation (3.3). The typical choice is taking as x a 
good approximation of a zero of / and as R a good approximation of (/'(x))“^, 
both obtained via a classic floating point algorithm, see for instance [9]. 

3.1. A Residual form for the Kraw^czyk operator. We now introduce 
the concepts that are needed to apply the modified Krawczyk method to solve 
a matrix equation such as (1.1). The Frechet derivative [16] of a Frechet 
differentiable matrix function F : ^ at a point X G ^ 

linear mapping Lp : such that for all i? G 

F{X + E)- F{X) - Lf{X,E) = o(||.F||). 

Since Lp is a linear operator, we can write 

Yec{Lp{X,E)) = Kp{X)vec{E), 

2 2 

for a matrix Kp(X) G C" that depends on L but not E. One refers to 
Kp{X) as the Kronecker form of the Frechet derivative of T at A. 

In the case of the continuous-time algebraic Riccati equation (1.1), we apply 
the Krawczyk method to the function E : ^ (^nxn (defined as 

E{X) := A*X + XA + Q- XGX, 
which appeared before in [14]. For this function, one has 
Lp(X, E) = E(A - GX) + {A* - XG)E. 

Lemma 2.1 part 5 turns out that its Kronecker form is 
Kp{X) = 4 (8) {A* - XG) + (A - GXf ® 4. 

When X = A*, we can write this expression in an alternate form as 

Kp{X) = In®{A- GX)* + {A- GX)'^ <S> In- (3.4) 

We wish to use the modified Krawczyk algorithm on the function obtained 
by regarding F as a vector map / : C'^, with N = n‘^, defined by 

/(x) := vec(A*A +AA + g- AGA), x = vec(A). (3.5) 

The following result, which is a slight variation of a theorem in [14], shows 
that the Frechet derivative can be used to obtain an enclosure for the slope in 
the modified Krawczyk method. We report it, with a different proof from the 
one in [14], because this presentation will be more convenient in the following 
development of our method. Due to this reformulation, we will get a weaker 
result with respect to uniqueness. 
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Theorem 3.3. Let X he an interval matrix, and Ki;’(X) = In® {A — GX)* + 
{A — GX)^ ® In he the interval evaluation of Kf{X) in (3.4). Then for each 
F, y' G X such that Y = Y*, it holds that S{f;y,y') G Kj7’(X), where y = 
vec{Y),y' = vec(y'). 

Proof. We have 

vec(F(y) - F{Y')) = vec((^* - yG)(y - Y') + (F - Y'){A - GY')) 

= vec((^ - Gy)*(y - Y') + (y - y'){a - gy')) 

= {In ® {A- GY)* + (A - GF')^ ® In) vec(y - y'), 
hence by the inclusion property of interval arithmetic 

S{f; y, y') = {In ® {A - GY)* + {A - GY')^ ® In) G Kf(X). (3.6) 

□ 

The next ingredient that we need to apply the Krawczyk algorithm is the 
matrix R. One would like to use R ~ {Kf{X))~^, where X = X* is an ap¬ 
proximation of the stabilizing solution to the CARE (1.1) computed in float¬ 
ing point arithmetic. However, this is the inverse of an x matrix, whose 
computation would cost 0{n^) floating point operations in general. Even con¬ 
sidering the Kronecker product structure of Kf{X), there is no algorithm in 

literature to compute R explicitly with less than 0{n^) arithmetic operations. 

2 

The action of R, that is, computing the product Rv given a vector v G C”' , can 
be computed with 0{n^) operations with methods such as the Bartels-Stewart 
algorithm [3]. However, this method cannot be used effectively in conjunction 
with interval arithmetic due to excessive wrapping effects, as argued in [10]. 

The work [14] (and, earlier, on a similar equation, [9]) contains an alterna¬ 
tive method to perform this computation with complexity 0{n^), in the case 
when A—GX is diagonalizable, where A is a numerical solution of CARE (1.1). 
Assume that an approximate eigendecomposition of A — GX is available, that 
is, 

A-GX ^ VAW with y, IF, A G C”""", (3.7a) 

A«diag(Ai,A 2 ,...,A„), VW^In. (3.7b) 

We write instead of = because F, IF « V~^ and A,, i = l,...,n are 
computed numerically with a standard method such as MATLAB’s eig. So, 
equality does not hold (in general) in the mathematical sense. Once these 
quantities are computed, we can factorize Kf{X) by replacing in the first 
term of (3.4) with and in the second term with W*InW~* and then 
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using Lemma 2.1, so that 

Kf{X) = In® {A- Gxy + {A- GXf ® In 

= {V-^ ® W*){In ® {W{A - GX)W-y* 

+ (V-^A - GX)Vf ® In){V^ ® W-*), 
and choose R as 

R = {V-^ ®W*)A-^{V^ ®W-*), (3.8) 

where A = 0 A* + ® R- Then, R ~ {Kf{X))~^ holds since if X is close 

enough to the stabilizing solution of (1.1), then one can expect that W(A — 
GX)W~^ and also V~^{A — GX)V to be close to A. So, the computation 
of an enclosure, vec(L), for I := —Rf{x) in Xf{x,R,z,S) can be done using 
exclusively the matrix-matrix operations, as shown in Lines 1-9 of Algorithm 1. 

For the latter term in each member of JCf{x,R,z,S), however, we get 

u := (42 - RS)z = (42 - {V-^ ® W*)A-^{V^ ® IF”*) 

(4 ® (A - GY f + {A- GY'f ® In))z 
= ®W*)A-^ 

{A -In® (WiA - GY)W-y* - {V-^A - GY')Vf ® 4) 
{V^ ®W-*))z, 

in which R 2 has replaced by ® W*W~*, and Y,Y' G X with Y = Y*. 

Then, Algorithm 2 Lines 2-6 will compute an enclosure for this term as the 
interval matrix U whose vectorization contains this term. 

Another point to note is that we can transform the multiplication T”^ vec(M), 
for an n X n matrix M and a diagonal matrix T, into M./N, where N is de¬ 
fined by Nij = Tjj-l-rjj, using point 6 of Lemma 2.1, and similarly for interval 
matrices using point 2 of Lemma 2.2. This point will appear in, for example. 
Algorithms 1 Line 8, and Algorithm 2 Line 6. 

The standard method [28] to obtain an interval vector z = vec(Z) that sat¬ 
isfies (3.3) is an iterative one. We start from the residual matrix Zq := F(X), 
that is, the interval evaluation of F{X), and proceed alternating successive 
steps of enlarging this interval with a technique known as e-inflation [28], ap¬ 
plying the Krawczyk operator to it, Zj+i = kj(x, i?, Zj, S). This procedure 
terminates when (and if) we find an interval for which (3.3) holds; it is ulti¬ 
mately a trial-and-error procedure, which is not guaranteed to succeed: the 
operator k j- may simply not contract its interval argument Zj sufficiently. This 
may be due to ill-conditioning of the original equation, to a bad choice of i?, or 
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to the excessive growth of the intervals in the numerical computations (wrap¬ 
ping effect). 

Several slightly different versions of the iterative procedure to obtain a valid 
interval for inclusion appear in literature; some involve intersecting the inter¬ 
vals obtained in different steps [9, 11, 15], and some involve two attempts at 
inclusion in each iteration [9, 15]. We use here the simplest approach, follow¬ 
ing [14, 28]. The exact strategy is shown in Algorithm 1 (and its subroutine 
Algorithm 2). The algorithm with these choices coincides with the algorithm 
presented in [14], except for the fact that [14] presents it for a generic Hermit- 
ian solution. 

In all algorithms, whenever the evaluation order of an expression is not 
specified exactly due to missing brackets, we evaluate from left to right. 


Algorithm 1 Computation of an interval matrix X containing a solution of 
CARE (1.1). 

I: Compute an approximate stabilizing solution X of CARE (1.1) using any 
floating point algorithm 

2: Compute approximations V, W, A for the eigendecomposition of A — GX 
in floating point {Eor instance, using the MATLAB command eig} 

3: Compute D := (Dij) such that Dij = An + Ajj 

4: Compute interval matrices ly and Ivu containing V~^ and W~^, respec¬ 
tively {For instance, using verifylss.m from INTLAB.} If this fails, or 
if D has any zero elements, return failure 
5: X = {X, 0) {To ensure that operations involving X are performed in a 
verified fashion with interval arithmetic} 

6: F = A*X -|- XA + Q — XGX {Using verified interval arithmetic} 

7: G = %FU 
8: H = G./D 
9: L = -lU*HIy 
10: Z = L 

II: for A: = 1, . . . ,kynax do 

12: Set Z = 0(0, Z • (1,0.1) -|- (0, realmin)) {e-inflation technique} 

13: Compute K using Algorithm 2 

14: if K C int(Z) {successful inclusion} then 

15: Return X = A -|- K 

16: end if 

17: Z = K 

18: end for 

19: Return failure {Maximum number of iterations reached} 
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Algorithm 2 Computation of an interval matrix K such that vec (K) = 
kf{x,R,z,S) encloses ICf{x,R,z,S). 

1: Input A, G, Q, X, Z 

{Additionally, in this subfunction we use V, W, ly, Iiy, A, D, L which are already 
computed in Algorithm 1} 

2: M = V^ZV 

3: N = W{A-G{X+ Z))lw 
4: 0 = Iy(A-G(X + Z))C 
5: P = (A - N)*M + M(A - O) 

6: Q = P./£I 
7: U = W*Qlv 
8: K = L + U 
9: Return K 


Notice the e-inflation, which is performed by adding (0, realmin) to the 
computed interval. Throughout the paper, realmin denotes the smallest pos¬ 
itive normalized floating point number. 

All the operations in Algorithm 1 are matrix-matrix computations requiring 
arithmetic operations, so its total cost is 0{n^s), where s is the number 
of steps needed before success. 


3.2. Affine transform enclosure. The main difficulty in using interval arith¬ 
metic to verify the existence of a solution is the so-called wrapping effect [28] : 
given interval matrices A and B, the set {AB : A G A,B G B} is not (in 
general) an interval matrix, and to represent it in interval form we have to 
enlarge it by replacing with an enclosing interval. The same effect happens 
with most interval operations, not only matrix multiplications, and it is more 
pronounced in presence of repeated successive operations, intervals with large 
radius, and ill-conditioned matrices. This increase may prevent us to verify 
computationally the critical condition (3.1). 

Reducing the impact of the wrapping effect can give a reduction in the 
diameter of the computed solution enclosure and also in the computational 
time, since it can reduce the number of iterations needed before a successful 
inclusion is computed. 

In this section we describe a technique for reducing the wrapping effect in 
the modified Krawczyk method, which has already been successfully applied 
to several matrix equations [9, 11]. The main idea is applying the verification 
algorithm to a modified function / obtained from / via an affine transforma¬ 
tion; in this way, we reduce the number of interval operations to perform inside 
the verification procedure. 
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Assuming that V and W defined in (3.7) are nonsingular, we define the 
function 

f(x) := (V^ 0 W-*)f{{V-^ (g) W*)x). (3.9) 

If X = vec(X) is an approximate solution to f(x) = 0, then x := (V'^iS>W~*)x 
is an approximate solution to f(x) = 0. The Kronecker form of its matrix 
formulation F{X) is given by 

Kp{X) = {V^ ® W-*)Kf{X){V-^ ® W*), X = W*XV-\ 

Moreover, let x = vec(X) := x + z, where z = vec(Z). A set of slopes for / on 
X can be dehned as 

5 := {S(f;y,y): y,y G x}. 

Dehning y := {V~'^ (g) W*)y, y' := {V~'^ (g) W*)y', we have 

S{f;y,y')iy - y') = Ry) - Ry') 

= iV^0W-*)if{y)-fiy')) 

= iV^0W-*)S{f-,y,y'){y-y') 

= {V^ ® W-*)S{f-, y, y'){V-^ ® W*){y - v')- 

Hence 

S{f-, y, y') = ® W-*)S{f- y, y'){V-^ ® IT*). 

In particular, if we combine this result with Theorem 3.3, we can take in the 
Krawczyk operator (3.2) 

S := S(/) = In®{W{A-GX)W-^y + {V-^{A-GX)vR®In. (3.10) 

where 

X = W*XV~\ X = l + Z, 

as long as X is Hermitian. 

Observe that 

In®iW{A-GX)W-y* + {V-\A-GX)Vf^In ~ 4(g)A* + A^(g)4, 
so a natural choice for R is the diagonal matrix 
R = A-\ 

in which A is defined as in (3.8). 

Now, we compute an enclosure for Xj{x,R,z,S) := {—RRx) + {In — 
RS)z,S £ S,z G z} which can be written as k.^{x, R,z,S) in which x is 

an approximate solution for (3.9), R is A“^, S = {S'(/;y,y'), y,y' £ x. := 
iy'^ (g) IT“*)x + z}, and z := vec(Z). As in Algorithm 1, we also take care 
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that the quantities which are not available exactly are enclosed into com¬ 
putable quantities in interval forms, for instance ly and I\y are interval matri¬ 
ces which are known to contain the exact value of V~^ and W~^, appropriately. 
More details for computing the superset 

kj(x, R, z, S) = —Rf{x) + (In — RS)z 

= -A-\{V^ 0W-*)f{x) 

- (A - 4 ® {W{A - 

for JCj{x,R,z,S), are displayed in Algorithm 4. The complete algorithm is 
shown in Algorithm 3. 


Algorithm 3 Computation of an interval matrix X containing a solution of 
CARE (1.1). 

1: Compute an approximate stabilizing solution X of CARE (1.1) using any 
floating point algorithm 

2: Compute approximations V, W, A for the eigendecomposition of A — GX 
in floating point {Eor instance, using the MATLAB command eig} 

3: Compute D := (Dij) such that Dij = An + Ajj 

4: Compute interval matrices ly and Iw containing V~^ and W~^, resp. 
{Eor instance, using verifylss .m from INTLAB.} If this fails, or if D has 
any zero elements, return failure 

5: X = {X, 0) {To ensure that operations involving X are performed in a 
verified fashion with interval arithmetic} 

6: F = A*X + Q + X(A - GX) 

7: F = I^FR 
8: L = -F./D 
9: Z = L 

10: for fc = 1 , . . . , krnax do 

11: Set Z = 0(0,^ • (1,0.1) -|- (0, realmin) {e-inflation technique} 

12: Compute K using Algorithm 4 (or Algorithm 5) 

13: if K C int(Z) {successful inclusion} then 

14: Return X = X W*Klv 

15: end if 

16: Z = K 

17: end for 

18: Return failure {Maximum number of iterations reached} 
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Algorithm 4 Evaluating K with vec(K) = z,S) encloses 

)C^{x, R, z, S). 

1: Input A, G, Q, X, Z 

{Additionally, in this sub-function we use V, W, ly, Iw, A, D, L which are already 
computed in Algorithm 3} 

2: M = W*Zlv 

3: N = I{^.(A-C?(A + M))*1E* 

4: 6 = Iy(A-G(X-LM))E 
5: P = (A* - N)M -f M(A - 6) 

6 : t = P./D 
7: K = L4-U 
8: Return K 


Note that computing L and U requires fewer dense n x n interval matrix 
multiplications than computing L and U, so the impact of the wrapping effect 
is reduced. This is the reason why one expects Algorithm 3 to work in more 
cases than Algorithm 1. 

An important observation is that the last transformation X = X -|- IT*KIy 
happens after the Krawczyk verification procedure. So, while the procedure 
guarantees that only one zero Xg of / is contained in W~*XV + K, when we 
return to the original setting and compute an enclosure for X = X -|- IT*KIy, 
other solutions of (1.1) may fall into this enclosure. Hence, Algorithm 3 alone 
does not guarantee that there is a unique solution of (1.1) in X, nor that this 
solution is the stabilizing one. We resolve with this issue in Section 3.5. 

Another small improvement introduced in this algorithm is gathering X in 
Line 6 of Algorithm 3, in order to reduce the wrapping effect. 


3.3. Verifying a different Riccati equation. Another possible modifica¬ 
tion to the verification process consists in modifying the equation into one 
with (possibly) better numerical properties. The idea stems from the formu¬ 
lation (1.2) of a CARE as an invariant subspace problem. We start from the 
following result. 


Lemma 3.4. [51 The stabilizinq solution X, of CARE (1.1) is the only matrix 
X E such that 


H 

'In 


'In 

pa 

II 

■ A -G' 


Xg 


Xg 

— ^ — 

-Q -A* 


E C 


2nx2n 


(3.11) 


for some Hurwitz stable matrix R. Moreover, it holds that R = A — GXg. 
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The subspace im 




is called stable invariant subspace of the matrix H. 


We use this formulation to relate the solution Xg to the one of a different 
CARE. The following result is a natural result of the literature on algebraic 
Riccati equations (see e.g. [5]), and the idea used here is certainly not original, 
but we prove it explicitly because we do not have a reference with this exact 
statement. 


Lemma 3.5. Let Xg be the stabilizing solution of (1.1). Suppose that P G 
£ 2 nx 2 n ^ nonsingular matrix such that P~^HP has the same structure as 
H, i.e., 


P-^HP 


Ap 

-Qp 


—Gp 

-A*p 


gC 


2nx2n 


(3.12) 


for some matrices Ap,Gp = G*p,Qp = Qp G Let Yg be the stabilizing 

solution of the CARE 


A*pY + YAp + Qp = YGpY, 

and Ui,U 2 G be defined by 

P W = K 

W ■ 

IfUi is invertible, then Xg = U 2 Uf~^. 
Proof. We have 


(3.13) 


P-^HP 


'In 


'In 





Rp, 


for the Hurwitz stable matrix Rp = Ap — GpYg. Multiplying both sides by P 
on the left we get 


H 


'Ui 


'u{ 



U2_ 


Rp, 


and then multiplying on the right by Ui 


-1 


H 


In 


In 

U2U{\ 


U2Ui\ 


UiRpU{ 


-1 


Since UiRpU^ ^ is Hurwitz stable, Lemma 3.4 gives us the thesis. 


□ 


The paper [23] contains a convenient strategy to construct a matrix P with 
a particularly simple form (a permutation matrix with some sign changes) for 






















VERIFIED CARE SOLUTIONS 


17 


which all the required assumptions hold and in addition Yg is bounded. Define 
for each k = 1,2,... ,n 


Sk := 


In Ekk 

Hkk 


E, 


kk 


In ~ E, 


kk 


G C' 


2nx2n 


where E^^ is the matrix which has 1 in position (k, k) and zeros elsewhere; in 
other words, swaps the entries k and n + A: of a vector in and changes 
sign to one of them. The matrices Sk are orthogonal and commute with each 
other. 


Theorem 3.6. [23, Theorem 3.4] Let X = { 11 , 12 ,... ,ik} he a subset of 

{1,2,..., n}, and P = Si,-- ■ Si^. Then 

(1) For each choice of I, the matrix P~^HP has the structure (3.12). 

(2) For each r > y/2, one can find I such that Ui is nonsingular and Yg has 
all its elements bounded in modulus by r (referring to the definitions 
of Ui and Yg in Lemma 3.5). 

These results suggest an alternative verification strategy: 

(1) Compute P satisfying Theorem 3.6. 

(2) Form the coefficients Ap,Gp and Qp, which can be obtained from the 
entries of H only using permutations and sign changes. 

(3) Using one of the various verification methods for CAREs, compute an 
interval Y containing the stabilizing solution Yg. 

(4) Compute 

In 

yJ ’ 

which, again, requires only rearranging the entries and changing their 
signs, and hence can be done without wrapping effects. 

(5) Compute using interval arithmetic a solution X to the linear system 

XUi = U 2 . 

Then, clearly, X contains the true stabilizing solution Xg of (1.1). Again, since 
the interval matrix X computed in the last step is only a solution enclosure 
and suffers from wrapping effect, it might be the case that other solutions of 
the CARE (1.1) are contained in X in addition to Xg. 

The MATLAB toolbox [27] contains algorithms to compute a subset X (and 
hence a matrix P) satisfying the conditions of Theorem 3.6, for every r > \/2, 
in time bounded by 0{n^ log.,, n). The factor log.,, n is a worst-case factor only, 
and in our experience for most matrices fine-tuning the choice of r does not 
have a big impact on neither performance nor stability. Here, we always use 
the method with its default value r = 3. 

With this method, one transforms the problem of verifying (1.1) into the 
one of verifying (3.13); this latter Riccati equation has a stabilizing solution 


Ui 

U2 


= P 
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Ys whose entries are bounded in modulus by r, hence one may expect that 
less cancellation can take place in the algorithms. While there is no formal 
guarantee that this must happen, in practice, in most cases the eigenvector 
matrix Vp oi Rp = Ap — GpYg has a lower condition number than the one 
V of A — GXg, as we report in the experiments (see Table 5 in the following), 
and verification of (3.13) is often easier than verification of (1.1). Ultimately, 
this is only a heuristic approach, though. 

Let us analyze the computational complexity of Algorithm 6. 

Theorem 3.7. Algorithm 6 requires at most 0{n^log^n + n^s) floating point 
operations, where s is the number of steps required by the inner verification 
algorithm in Line 5. 

Proof. Computing X in Line 2 requires 0{n^) operations, using for instance 
the algorithm mentioned in [25] (based on the ordered Schur form of H and an 
additional Newton step with the residual computation performed in emulated 
quadruple-precision arithmetic). Forming P in Theorem 3.6 via the approach 
explored in [27] costs 0(n^log.,-n) floating point operations. Computing Y 
by using Algorithm 3 has cost 0{rfi) per step (and the same will hold for 
Algorithm 8 that we will introduce later): the cost for the eigendecomposition 
and the enclosures ly and Iw is again cubic in n, and all the other matrix- 
matrix operations (including the Hadamard divisions) in Algorithms 3 and 5 
have again cost 0{n^) at most, as they only involve n x n matrices. □ 

3.4. A new superset. According to Theorem 3.2, the computed interval 
matrix is guaranteed to contain a unique solution if the set S contains the 
slopes 5(/; y, y') for all y, y' G x. On the other hand, if we employ an interval 
matrix containing only the slopes S{f]x,y') for all y' G x, existence can be 
proved, but not uniqueness. Since we have already decided to forgo (for now) 
uniqueness, it makes sense to let go of it also when choosing the superset S. 

A simple modification to our proof of Lemma 3.3 gives a tighter inclusion 
for the slope superset by reducing the wrapping effect. 

Theorem 3.8. Let f be as in (3.5), X G be an interval matrix, and 

X gX be Hermitian. Then, the interval matrix 

/„ (8) (A - GXy + {A- GXf ® Ln 

contains the slopes S{f,x,y') for each y' G X where x = vec(X) and y' = 
vec(Y'). 


Proof. We repeat the proof of Lemma 3.3, with y = x, and replace the term 
Kj 7’(X) in (3.6) with the tighter inclusion (I„ (g) (A — GX)* -|- (A — GX)^ (g) 
In). □ 
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As a consequence of Theorem 3.8, we can replace (3.10) with 

S = /„ (8) (WiA - GX)W-^)* + {V-\A - GX)Vf 0 4 (3.14) 

in our modified Krawczyk algorithm applied to /, and it will still yield an 
interval matrix containing a (possibly non-unique) solution of (1.1). 


Algorithm 5 Evaluating K with vec(K) = k^(x,i?, z,S) encloses 

ICj:{x, R, z, S) with a tighter superset that does not guarantee solution unique- 
ness. _ 

{This algorithm is identical to Algorithm 4, apart from Line 3 which is 
replaced by the following} 

3: N = ryy{A-GX)*W* 


Algorithm 6 Computation of an interval matrix X containing a solution 
of (1.1) using permuted Riccati bases. 


I: 

2 : 

3: 

4: 

5: 

6 : 

7: 

8 : 


Input A, G, Q 

Compute an approximate stabilizing solution X of (1.1) using any floating 
point algorithm 

Compute a matrix P satisfying point 2 of Theorem 3.6 {For instance with 
the toolbox [27]} 

Compute Ap,Gp,Qp satisfying (3.12) 

Compute a verified solution Y to (3.13) using either Algorithm 3 or Algo¬ 
rithm 8 . If the verification fails, return failure 


Set 


Ui 

U2 


= p 


Compute an enclosure X for the solution of the interval system XUi = 
U 2 {For instance, using verifylss from INTLAB}. If this fails, return 
failure 
Return X 


3.5. Verification of uniqueness and stabilizability. As noted before, the 
modifications to the Krawczyk method introduced here do not ensure that 
the found interval matrix contains only one solution to (1.1). However, the 
following result holds. 

Theorem 3.9 ([6, Theorem 23.3]). The CARE (1.1) has at most one stabi¬ 
lizing solution. 
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A proof using the facts in [5] can be obtained by considering the eigenvalues 
of H. Let Xg be a stabilizing solution, and let Ai, A 2 , • • •, be the eigenvalues 
of A — GXg. Because of the formula (1.2), Ai, A 2 , • • •, A„ are also eigenvalues 
of H (see [5, Section 2.1.1]). Moreover, the eigenvalues of H have Hamiltonian 
symmetry, (see [5, Section 1.5]), so there are n more eigenvalues with positive 
real part. We have identified 2n eigenvalues of H, counted with multiplicity, 
and none of them is purely imaginary; hence H has no purely imaginary 
eigenvalue and [5, Theorem 2.17] holds. 

Hence, if all matrices contained in A — GX are Hurwitz stable, then the 
solution A G X is verified to be stabilizing (hence Hermitian) and unique. To 
verify stability, we can use the method described in [24], which is summarized 
in [25, Lemma 2.4]. The resulting method is described in Algorithm 7. In 
the algorithm, we use the notation to mean the real part of the complex 
number z. 


Algorithm 7 Verifying the Hurwitz stability of an interval matrix M. 

1: Input M 

2: Compute approximations V, W, A for the eigendecomposition of mid(M) in float¬ 
ing point {For instance, using the MATLAB command eig} 

3: V = (V, 0) 

4: i? = mag(lT(MV - VA)) 

5: S'= mag(/„ - VIT) 

6: e = the n x 1 matrix with = 1 for each i 
7: u = Re 
8: t = Se 

9: p, = max(u./(e — t)) 

10: r = u + fit 

11: if (max(t) < 1 and r + max(5l7(diag(A)))e < 0) then 

12: Return success {Every matrix M G M is Hurwitz stable} 

13: else 

14: Return failure 

15: end if 


Notice one subtle point: when we apply Algorithm 7 to A — GX we recom¬ 
pute V, A and W from the eigendecomposition of mid(A — GX); this differs 
slightly from using the values computed previously, which came from the eigen¬ 
decomposition of A — GX (because X was not available at that point). This 
choice gives better results in our experiments. The cost for this verification is 
again 0{n^) floating point operations. 

Hence, if the verification in Algorithm 7 succeeds for the solution enclosure 
X returned by either Algorithm 1 or Algorithm 3, then X contains exactly 
one solution of (1.1), and it is the stabilizing one. 
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4. A DIRECT FIXED-POINT METHOD 

While the methods described in the previous sections work for many exam¬ 
ples of Riccati equations, an essential limitation is that all of them require the 
closed-loop matrix A—GX to be diagonalizable. Products with the eigenvector 
matrix V and its inverse are required along the algorithm, and if these are ill- 
conditioned then the wrapping effects are more pronounced and the required 
inclusion K C int(Z) or K C int(Z) is less likely to hold. A striking example 
of this phenomenon is the first example in the benchmark set [4]. This is a 
simple 2x2 problem which appears in [4] as nothing more than a “warm-up 
example”, and yet all the verification methods described here (including those 
from [14] and [25]) fail. 

To solve this issue, we would like to propose a different method for verifi¬ 
cation. The procedure is based on some ideas which appear in the context of 
ADI methods [30]. While this method is somehow more primitive and works 
on a lower number of examples, it does not require that the closed-loop matrix 
be diagonalizable. 

We rewrite the CARE (1.1) as follows. Given any Hermitian X € one 

can write the exact stabilizing solution Xg of the CARE (1.1) as Xg = X -h Zg 
for an unknown Hermitian correction matrix Zg, and rewrite (1.1) as a Riccati 
equation in Z, 

A*Z+ZA+Q = ZGZ, with A = A-GX, Q = A*X+XA+Q-XGX. 

(4.1) 

The stabilizing solution of this equation is Zg, since A — GZg = A — GXg is 
Hurwitz stable. Note that the degree-two coefficient G is unchanged. Eor any 
s G C such that A — sin is nonsingular, (4.1) is equivalent to the fixed point 
equation 

Z = (A- sIn)~*{ZGZ -Q-Z{A + sin)). 

Thus, if we find an interval Z such that {A —sIn)~*{'ZG7j — Q— 'L{A +sin)) C 
Z, it follows from the Brouwer fixed-point theorem that (4.1) has a solution 
Z^ G Z, and that (1.1) has a solution A* G A -|- Z. 

This simple iterative method is effective when [A — sIn)~*'Zi[A + sin) does 
not suffer excessively from wrapping effects, since we can expect Q and the 
quadratic term ZGZ to be small. 

Are there any preconditioning transformations that we can make to improve 
the method? A possibility is applying a change of basis to the whole problem. 
Let V G be invertible; we set 

Zy = Av = V-^AV, Qv = V*QV, Gy = R-^GE"*, (4.2) 
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SO that (4.1) is transformed into 

AyZv + ZyAy + Qy = ZyGyZy. 

Continuing as above, we obtain the fixed-point equation 

Zy = {Ay — sIn)~*{ZyGyZy — Qy — Zy{Ay + sIn))- (4.3) 

If A is diagonalizable, we can set V as its computed approximate eigenvector 
matrix, as in (3.7). One can see then that the resulting method has several 
steps in common with the Krawczyk method described in the previous sections. 
This time, though, we are free to choose the matrix V without the risk of our 
method turning into a O(n^) one, since everything in (4.3) is computable 
explicitly with standard linear algebra operations. 

Some heuristic experimentation led us to the following choices: we take s 
equal to the approximation of — min{3ftA : A is an eigenvalue of A} computed 
in floating-point arithmetic (motivated by the idea to make Ay + sin small 
and Ay — sin large), and V as the orthogonal factor of the (computed) Schur 
factorization of ^4 ~ VTV~^ (motivated by the idea to concentrate most of 
the “weight” of V~^AV on its diagonal). The matrix A is an approximation of 
A — GXg, which is Hurwitz stable, so in exact arithmetic we would have s > 0 
and Ay —sin = V~^{A—sIn)V invertible, since its eigenvalues are Xi—s, where 
Aj are the eigenvalues of ^ — GXg, and thus have strictly negative real part. 
Hence these properties are likely to hold also for its computed approximation 

i. 

The resulting algorithm is described in Algorithm 8. 

Theorem 4.1. Algorithm 8 has a cost ofO{n^s) arithmetic operations, if the 
verification succeeds in s steps. 

Proof. Again, all the required operations in every step are matrix-matrix op¬ 
erations between n x n matrices. The Schur decomposition requires 0{n^) 
operations as well, in practice [12]. □ 

Once again, uniqueness is not guaranteed, but it can be deduced a posteriori 
if the verihcation of the stabilizing property of the computed inclusion interval 
X succeeds. 


5. Numerical experiments 

This section presents numerical experiments to validate the algorithms. We 
compare four different approaches: 

(1) The modified Krawczyk approach described in [14] and in Section 3.1. 
This corresponds to Algorithm 1. When the algorithm is successful, we 
check afterwards whether A — GX is Hurwitz stable using Algorithm 7. 
We call this approach Method H in the following. 
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Algorithm 8 Computation of an interval matrix X containing a solution 
of (1.1) using a simple fixed-point algorithm. 

I: Input A, G, Q 

2: Compute an approximate stabilizing solution X of (1.1) using any floating 
point algorithm 

3: Compute A (in floating point) as in (4.1) 

4: Choose s and V, for instance s ~ — min{5RA : A is an eigenvalue of A} and 

V as the (approximate) orthogonal Schur factor of A 

5: Compute an interval matrix ly containing V~^ {For instance, using 
verifylss from INTLAB} 

6: Compute interval matrices Ay, Gy, Qy containing Ay, Gy, Qy, respec¬ 
tively {Replacing X and V in (4.1) and (4.2) with X = {X,0) and 

V = (F, 0), respectively} 

7: Compute an interval matrix containing (Ay — sln)~^ {For instance, 
using verifylss from INTLAB} 

8: Set A; = 0 and Zy = —RQy 
9: for /c = 1, . . . ,kmax do 

10: Set Zy = 0(0, Zy • (1,0.1) -|- (0, realmin)) {e—inflation technique} 

11: Set Y = Is( — Qy — Zy(Ay -|- sIn — GyZy)) 

12: if Y C int(Zy) then 

13: Return X = A -|- lyZyly 

14: end if 

15: end for 

16: Return failure {Maximum number of iterations reached} 


(2) The method described in [25] (using the MATLAB implementation 
Mn.m published by its author). The method already includes running 
Algorithm 7 to check if the computed solution is Hurwitz stable, so we 
do not need any additional steps. We call this procedure Method M. 

(3) Algorithm 6, choosing as its subroutine to solve the transformed CARE (3.13) 
the Krawczyk-based Algorithm 3 and the modified superset trick used 

in Algorithm 5. This is a combination of all the improvements to 
Method H described in Section 3. We call this procedure Method K 
(where K stands for Krawczyk). When the algorithm is successful, we 
check afterwards whether A — GX is Hurwitz stable using Algorithm 7. 

(4) Algorithm 6 again, but using the fixed-point Algorithm 8 to solve the 
transformed CARE (3.13). This is a combination of the techniques 
described in Sections 3.3 and 4. We call this procedure Method F 
(where E stands for fixed-point). When the algorithm is successful, we 
check afterwards whether A — GX is Hurwitz stable using Algorithm 7. 
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The algorithms were tested in MATLAB 2015b with INTLAB v6, using unit 
round off tt = 2“®^ ~ 1.1 x 10“^®, and run on a computer with an Intel core 
Duo 2.66GHz processor and 6GB main memory. 

The required stabilizing solutions of GAREs are computed using the method 
described in [25] (ordered Schur method followed by one step of Newton re¬ 
finement in simulated quadruple precision). 

In order to assess the quality of the enclosures computed in each experiment 
we use the norm-wise relative error nre and the geometric average relative 
precision garp. The first error measure is dehned as 


nre(X) := mag 


|rad(X)|| 

IIXIL 


F 


This is the simplest possible bound for the (norm-wise) relative error 


||A,-mid(X)||^ 

obtained by taking mid(X) as an approximation of the solution. 

Following previous works [10, e.g.], we also report a component-wise error 
indicator garp based on the relative precision of an interval, rp(x), defined as 


rp(x) := min(relerr(x), 1), 


where relerr is the relative error of the interval x 
by 


relerr (x) 


rad(x) 
mid(x) ’ 

rad(x). 


if 0 ^ X, 
otherwise. 


(mid(x),rad(x)) defined 


We dehne our residual measure as the geometric average of rp(Xjj) 
garp(X) := 

The quantity — log (rp(x)) can be interpreted as the number of known correct 
digits of an exact value contained in x; so, loosely speaking, — log (garp(X)) 
represents the average number of known correct digits [11]. 


n w(Xij) 




X = (Xi,). 


5.1. CAREX Benchmark problems. We ran these algorithms on all the 
equations from the benchmark set described in [7], which contains experiments 
taken from the test suite GAREX [4], run with both default and non-default 
arguments. The results are reported in Tables 1-4, and a visualization of the 
results is in Figure 6. 

The Experiment number follows the order used in [7]. Note that this set 
of problems is designed to be challenging for non-verified CARE solvers in 
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machine arithmetic, so it is not surprising that the verification algorithms 
cannot deal with all of them with perfect accuracy. 

When the algorithms are successful, we report in Tables 1-3 the number 
k of required iterations of the outer Krawczyk loop. If the algorithm breaks 
down or does not converge within the maximum number of steps (which we 
took to be 50 for the iterative Methods H, K and F), then we write a star in 
the corresponding column. Method M is not iterative, therefore for it we put 
— in the column containing the number of iterations. 

The size of the problem (value of n) and the total time (in seconds, including 
the time required to verify the stabilizing property) taken on our test machine 
are reported, too, as well as the norm-2 condition number of V (used by 
Methods H, M and K) and the same quantity for the eigenvector matrix Vp 
of the closed-loop matrix Ap — GpYg used in the two Methods K and F. All 
these details are given in Tables 1-3 in the column named Problem property. 


Table 1. Comparison between various proposed methods 


Experiment 

Problem 

property 

Method H 

Method M 

Method K 

Method F 

number 

size 

nre 

k 

nre 

k 

nre 

k 

nre 

k 

in [7] 

cond(T) 

cond(Vp) 

garp 

time 

garp 

time 

garp 

time 

garp 

time 

1 

2 

NaN 

* 

NaN 

- 

NaN 


3.75e-15 

2 


7.75e+07 

3.17e+07 

NaN 

% 

NaN 

* 

NaN 

* 

4.18e-15 

2.93e-02 

2 

2 

9.67e-14 

1 

4.65e-15 

- 

1.21e-14 

1 

l.OOe-14 

3 


l.Ole+01 

1.15e+00 

1.04e-13 

2.00e-02 

4.97e-15 

7.59e-03 

1.27e-14 

2.28e-02 

1.06e-14 

2.70e-02 

3 

4 

3.93e-14 

1 

2.99e-15 

- 

3.70e-14 

1 

8.04e-14 

6 


9.73e+00 

5.11e+00 

2.80e-14 

2.28e-02 

2.12e-15 

1.03e-02 

5.02e-14 

2.90e-02 

1.05e-13 

4.28e-02 

4 

8 

1.02e-14 

1 

2.34e-15 

- 

7.76e-14 

1 

1.03e-13 

15 


1.23e+00 

2.18e+00 

1.49e-14 

1.74e-02 

3.42e-15 

9.05e-03 

1.05e-13 

2.44e-02 

1.38e-13 

5.72e-02 

5 

9 

6.73e-14 

1 

l.lOe-14 

- 

4.34e-13 

1 

2.06e-12 

42 


7.54e+01 

6.52e+01 

4.34e-14 

1.78e-02 

1.05e-14 

9.37e-03 

7.57e-13 

2.45e-02 

4.70e-12 

1.23e-01 

6 

30 

4.79e-13 

2 

3.35e-14 

- 

9.20e-09 

2 

NaN 

* 


l.lle+05 

3.48e+03 

2.92e-ll 

5.15e-02 

1.87e-12 

1.88e-02 

1.15e-08 

6.64e-02 

NaN 

* 

7 

2 

2.35e-16 

2 

2.12e-16 

- 

5.57e-16 

1 

7.47e-16 

3 


1.62e+00 

3.31e+00 

5.48e-16 

2.18e-02 

4.36e-16 

7.58e-03 

6.17e-16 

2.27e-02 

8.72e-16 

2.69e-02 

8 

2 

3.22e-08 

1 

1.78e-08 

- 

3.67e-16 

1 

8.55e-16 

2 


l.Ole+00 

2.42e+00 

4.75e-10 

1.68e-02 

3.04e-10 

1.08e-02 

4.83e-16 

2.34e-02 

1.88e-10 

2.52e-02 

9 

2 

6.41e-16 

1 

1.92e-16 

- 

3.34e-10 

1 

2.25e-09 

7 


1.22e+00 

6.99e+01 

2.59e-15 

1.69e-02 

4.87e-16 

8.31e-03 

3.36e-10 

2.28e-02 

2.26e-09 

3.64e-02 

10 

2 

NaN 

% 

5.36e-12 

- 

3.00e-08 

1 

1.59e-08 

45 


6.80e+01 

l.lle+00 

NaN 

% 

5.36e-12 

8.49e-03 

3.00e-08 

2.15e-02 

1.59e-08 

1.34e-01 

11 

2 

9.65e-16 

1 

6.29e-16 

- 

2.45e-15 

1 

4.53e-15 

2 


3.74e+00 

l.Ole+00 

1.04e-15 

2.20e-02 

6.75e-16 

7.71e-03 

2.63e-15 

2.90e-02 

4.87e-15 

3.00e-02 
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Table 2. Comparison between various proposed methods 


Experiment 

Problem 

property 

Method H 

Method M 

Method K 

Method F 

number 

size 

nre 

k 

nre 

k 

nre 

k 

nre 

k 

in [7] 

cond(T) 

cond(Vp) 

garp 

time 

garp 

time 

garp 

time 

garp 

time 

12 

2 

7.96e-16 

1 

3.22e-16 

- 

1.53e-15 

1 

3.90e-ll 

2 


1.42e+03 

l.Ole+00 

8.94e-16 

2.20e-02 

3.74e-16 

7.50e-03 

2.29e-15 

2.82e-02 

6.56e-ll 

2.95e-02 

13 

2 

9.41e-09 

1 

2.23e-09 

- 

6.99e-16 

1 

NaN 

* 


2.42e+00 

l.Ole+00 

3.01e-10 

1.69e-02 

7.20e-ll 

8.29e-03 

1.15e-15 

2.28e-02 

NaN 

* 

14 

2 

1.65e-15 

1 

2.69e-16 

- 

9.49e-16 

1 

4.37e-15 

3 


l.Ole+00 

l.Ole+00 

1.92e-15 

1.65e-02 

3.14e-16 

7.24e-03 

1.09e-15 

2.29e-02 

5.09e-15 

2.71e-02 

15 

2 

4.72e-ll 

1 

3.36e-12 

- 

2.92e-15 

1 

NaN 

* 


l.Ole+00 

1.29e+00 

4.72e-ll 

1.67e-02 

3.36e-12 

8.09e-03 

2.39e-15 

2.27e-02 

NaN 

* 

16 

2 

NaN 

% 

5.87e-10 

- 

9.38e-12 

1 

NaN 

* 


l.OOe+00 

1.29e+00 

NaN 

% 

5.87e-10 

8.07e-03 

5.63e-12 

2.15e-02 

NaN 

* 

17 

2 

2.42e-15 

1 

2.23e-16 

- 

4.80e-15 

1 

1.25e-14 

5 


l.Ole+00 

2.62e+00 

2.69e-15 

2.23e-02 

2.23e-16 

9.29e-03 

4.98e-15 

2.83e-02 

1.35e-14 

3.92e-02 

18 

2 

NaN 

* 

NaN 

- 

NaN 


NaN 



l.Ole+00 

2.62e+00 

NaN 

% 

NaN 

* 

NaN 

* 

NaN 

* 

19 

3 

3.53e-15 

1 

2.73e-16 

- 

3.67e-15 

1 

6.64e-15 

3 


l.Ole+00 

l.Ole+00 

1.20e-14 

1.98e-02 

8.76e-16 

7.88e-03 

1.19e-14 

2.29e-02 

2.16e-14 

2.72e-02 

20 

3 

9.95e-05 

3 

3.87e-05 

- 

4.77e-15 

1 

1.73e-14 

3 


l.Ole+00 

l.OOe+00 

1.18e-04 

2.78e-02 

4.53e-05 

8.36e-03 

6.08e-12 

2.30e-02 

2.08e-14 

2.70e-02 

21 

4 

1.29e-14 

1 

4.76e-15 

- 

1.26e-13 

1 

3.81e-13 

11 


9.01e+00 

3.58e+00 

1.73e-14 

2.25e-02 

6.41e-15 

8.79e-03 

1.39e-13 

2.92e-02 

4.32e-13 

5.84e-02 

22 

4 

7.28e-05 

2 

3.52e-06 

- 

1.70e-04 

1 

NaN 



1.22e+01 

5.94e+00 

4.81e-06 

3.02e-02 

2.96e-07 

1.04e-02 

8.34e-06 

2.90e-02 

NaN 

* 
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Table 3. Comparison between various proposed methods 


Experiment 

Problem 

property 

Method H 

Method M 

Method K 

Method F 

number 

size 

nre 

k 

nre 

k 

nre 

k 

nre 

k 

in [7] 

cond(T) 

cond(Vp) 

garp 

time 

garp 

time 

garp 

time 

garp 

time 

23 

4 

3.25e-14 

1 

1.78e-15 

- 

3.43e-14 

1 

1.30e-13 

11 


1.43e+01 

1.77e+00 

3.00e-14 

2.23e-02 

1.71e-15 

7.96e-03 

4.61e-14 

2.92e-02 

1.73e-13 

5.84e-02 

24 

4 

NaN 

* 

NaN 

- 

NaN 

* 

NaN 

* 


1.74e+00 

1.74e+00 

NaN 

* 

NaN 

* 

NaN 

* 

NaN 

* 

25 

77 

4.08e-12 

1 

3.66e-13 

- 

4.70e-ll 

1 

3.16e-10 

12 


4.98e+01 

1.87e+01 

3.50e-ll 

2.11e-01 

3.14e-12 

1.13e-01 

2.64e-10 

3.03e-01 

1.63e-09 

4.34e-01 

26 

237 

4.75e-ll 

1 

4.35e-12 

- 

2.27e-09 

1 

1.26e-08 

17 


2.41e+02 

8.93e+01 

8.89e-10 

3.05e+00 

8.21e-ll 

1.53e+00 

1.41e-08 

4.56e+00 

7.09e-08 

6.79e+00 

27 

397 

NaN 

* 

6.71e-12 

- 

8.73e-09 

1 

6.69e-08 

19 


1.31e+02 

4.84e+01 

NaN 

* 

2.27e-10 

5.26e+00 

5.50e-08 

1.82e+01 

3.87e-07 

2.87e+01 

28 

8 

4.35e-15 

1 

1.67e-15 

- 

4.35e-15 

1 

1.30e-14 

4 


l.Ole+00 

l.Ole+00 

8.95e-15 

1.74e-02 

3.44e-15 

7.72e-03 

8.95e-15 

2.37e-02 

2.66e-14 

3.04e-02 

29 

64 

4.12e-13 

1 

4.99e-14 

- 

4.12e-13 

1 

7.12e-13 

4 


l.Ole+00 

l.Ole+00 

1.97e-07 

4.56e-02 

2.38e-08 

3.63e-02 

1.97e-07 

7.01e-02 

3.40e-07 

8.96e-02 

30 

21 

NaN 

* 

NaN 

- 

3.90e-04 

1 

NaN 

38 


2.42e+09 

2.78e+00 

NaN 

* 

NaN 

* 

3.79e-04 

4.34e-02 

NaN 

2.08e-01 

31 

21 

NaN 

* 

NaN 

- 

NaN 

1 

NaN 

* 


2.42e+09 

2.88e-l-02 

NaN 

* 

NaN 

* 

NaN 

5.06e-02 

NaN 

* 

32 

100 

6.57e-12 

1 

1.13e-12 

- 

6.57e-12 

1 

NaN 

* 


l.Ole+00 

l.Ole+OO 

2.27e-ll 

1.52e-01 

3.90e-12 

1.36e-01 

2.27e-ll 

2.33e-01 

NaN 

* 

33 

60 

2.04e-14 

1 

2.67e-13 

- 

2.77e-10 

1 

NaN 

* 


1.91e+01 

1.55e+01 

4.67e-14 

1.12e-01 

6.10e-13 

4.97e-02 

4.17e-10 

1.58e-01 

NaN 

* 
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Figure 1. nre of Method H vs. cond(F) 



Figure 2. nre of Method M vs. cond(yp) 
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Figure 3. nre of Method K vs. cond(F) 
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Figure 4. nre of Method F vs. cond(F) 




VERIFIED CARE SOLUTIONS 


31 



Figure 5. cond(Vp) vs. condCF). Most of the points lie be¬ 
low the axes bisector (drawn in red), which means that the 
condition number of Vp is generally lower than the one of V. 
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Figure 6. Values of garp for each experiment number. A full bar means that the method failed 
to compute an enclosure. 
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Table 4 reports the result of checking the stabilization property; a plus sign 
means that the property is verified, a minus sign means failure to verify the 
property, and a star means that the algorithm had already failed to compute 
an inclusion interval. As one can see, there is only a very limited number of 
cases in which the stabilization procedure fails. 

Remarks are in order on some of the problems. 

Experiment 1: This is an example of the phenomenon described in the 
beginning of Section 4: the closed-loop matrix A — GXg is not diago- 
nalizable. The coefficient matrices for this example are 


A = 


0 1 
0 0 




0 0 
0 1 


and Q 


1 0 
0 2 


The exact value of the closed loop matrix for the original and trans¬ 
formed equations are respectively 


A-GX, 


0 1 
-1 -2 


and Ap — GpYg 


-2/3 1/3 ■ 

-1/3 - 4 / 3 J ’ 


both with a double (defective) eigenvalue in —1. The approximation 
X computed with the Schur method satisfies ||Xs — X\\ = 1.68e — 15. 
The matrix A — GX is diagonalizable with two very close eigenvalues. 
Hence, the computed condition numbers of V and Vp are both large, 
and the first three algorithms, which are based on the diagonalization 
of an approximation of A — GXg, fail. On the other hand, the fixed- 
point algorithm does not encounter any difficulty and returns a tight 
interval X containing the stabilizing solution. The condition number 
of the eigenvector matrix of mid(A — GX) is 7.75e7, but the verification 
with Algorithm 7 succeeds nevertheless. 

Experiments 30 and 31: In Method F for problem 30 and Method K 
for problem 31, we report termination in a finite number of itera¬ 
tions, but NaN for the error. In these problems, the verification al¬ 
gorithm succeeds for the Riccati equation (3.13), but the resulting in¬ 
terval Y cannot be converted into a solution interval X for (1.1) using 
Lemma 3.5, because the interval matrix Ui computed as described in 
Section 3.3 contains singular matrices, hence the solution set X is un¬ 
bounded. So the method fails to produce a solution enclosure for (1.1). 

Another interesting observation is that Method K needs only one iteration 
in all experiments when it works apart from one case (Experiment 6), i.e., the 
crucial relation (3.3) is already fulfilled for fe = 1 in all the other cases. 

When they are successful. Methods H and K are comparable with respect to 
execution time as well as with respect to the quality of the enclosure. However, 
there are cases in which Method H is not successful, and this comprises cases 
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Table 4. Results for stabilizing property in all methods 


Experiment number 

Method H 

Method M 

Method K 

Method F 

1 


* 

* 

+ 

2 

+ 

+ 

+ 

+ 

3 

+ 

+ 

+ 

+ 

4 

+ 

+ 

+ 

+ 

5 

+ 

+ 

+ 

+ 

6 

+ 

+ 

+ 


7 

+ 

+ 

+ 

+ 

8 

+ 

+ 

+ 

+ 

9 

+ 

+ 

+ 

+ 

10 


+ 

- 

- 

11 

+ 

+ 

+ 

+ 

12 

+ 

+ 

+ 

+ 

13 

+ 

+ 

+ 


14 

+ 

+ 

+ 

+ 

15 

+ 

+ 

+ 

* 

16 


+ 

+ 


17 

+ 

+ 

+ 

+ 

18 





19 

+ 

+ 

+ 

+ 

20 

+ 

+ 

+ 

+ 

21 

+ 

+ 

+ 

+ 

22 

+ 

- 

+ 

* 

23 

+ 

+ 

+ 

+ 

24 

* 


* 


25 

+ 

+ 

+ 

+ 

26 

+ 

+ 

+ 

+ 

27 


+ 

+ 

+ 

28 

+ 

+ 

+ 

+ 

29 

+ 

+ 

+ 

+ 

30 

* 

* 

- 

- 

31 



- 


32 

+ 

+ 

+ 


33 

+ 

+ 

+ 



































VERIFIED CARE SOLUTIONS 


35 


with small dimensions (e.g. 2 in Example 10) as well as cases with large 

dimensions (e.g. 397 in Example 27). 

Method M is significantly faster than the other algorithms. We remark, 
though, that MATLAB, being an interpreted language, is often not reliable 
in evaluating computational times. In particular, INTLAB is implemented 
entirely in MATLAB code, and its running time does not always match the 
theoretical complexity, especially when dealing with small matrices. Eor Meth¬ 
ods K and E, which rely on Algorithm 6, another consideration is that the 
computation of the matrix P using the toolbox [27] requires in its default 
implementation a tight double for loop on the matrix entries. MATLAB ex¬ 
ecutes loops of this kind much more slowly than operations on full matrices; 
hence comparing running times may show a larger discrepancy than the actual 
difference in performance between the algorithms. 

Methods K and M are the most reliable, and fail only on very ill-conditioned 
examples. Interestingly, the errors obtained by the two approaches differ by 
orders of magnitude on several problems, in both directions; there are also 
examples in which either one fails while the other succeeds. So there is no 
clear winner among the two. 

Method E has the largest number of failures. Despite that, it is useful in 
special cases (such as in Experiment 1) in which the other algorithms have 
difficulties, particularly when the closed-loop matrix is not diagonalizable. 

In many of the examples the performance of the methods based on diagonal¬ 
izing the closed-loop matrix is (loosely) related to the condition number of V 
(or Vp, when it is used). To visualize this relationship, we show in Eigures 1-4 
scatter plots of the obtained accuracy vs. the value of this condition number in 
the various examples. When the magnitude of cond(E) is moderate, cond(I/p) 
has typically the same order of magnitude, but in some cases when cond(E) is 
large cond(yp) seems to be considerably lower, as shown in Eigure 5. There is 
only one case in which cond(yp) is considerably larger than cond(E), that is. 
Experiment 9 (1.22 vs. 69.8). This shows experimentally that switching from 
the formulation (1.1) to (3.13) is benehcial. 

5.2. Experiments with varying sizes. In view of the fact that the three 
Methods H, K and E are iterative taking an unspecified number of steps, and 
that the last two require a factorization which may require 0{n^ log.,- n) in the 
worst case, when n is the size of X in (1.1), the reader may wonder how the 
time taken by the various algorithm scales with the dimension n in practice. 
We have tested all algorithms on [4, Problem 15], which is a problem designed 
explicitly to check how Riccati solvers scale with the dimension of the equation. 
We have generated the test problem in 30 different sizes equally distributed in 
logarithmic scale between 10 and 1000, and we have tested the four algorithms 
on these examples. The resulting CPU times are reported in Eigure 7. 
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n 


Figure 7. CPU times for verification on a scaled version of 
Experiment 15 vs. dimension n. 

Overall, the results shows that all methods scale essentially with 0{n^), and 
in particular that Methods K and F stay within a moderate factor of the time 
taken by Method M. In the two largest experiments n = 853, n = 1000, Method 
K is the only one to succeed: Method M fails, while Method F delivers a 
solution enclosure for which the stabilizing property cannot be proved. Method 
H fails for each n > 204. Verification of the stabilizing property succeeds in 
all other cases apart from the two mentioned above for Method F. 

The MATLAB code used for the experiments is available online on https: //bitbucket. org/f ph/ 

6. Summary and Outlook 

We have introduced several improvements to the method in [14], borrowing 
ideas from both the interval arithmetic and the matrix equations literature. 

The resulting method has been tested on several standard benchmark experi¬ 
ments, and is competitive with the one introduced in [25], returning a smaller 
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solution enclosure in several of the experiments. Moreover, the new fixed-point 
method described in Section 4 is a useful addition to the battery of existing 
verification methods; it is especially useful in the cases in which the closed-loop 
matrix is not diagonalizable. 

There is no single algorithm that beats all the others on all the benchmark 
problems; hence it is important to have several methods available, each with 
its strengths and drawbacks. Overall, all but two of the problems in this 
challenging set of experiments could be verified with success. 

A number of open problems remain: first of all reducing to zero the num¬ 
ber of remaining failures in the methods. Of particular interest would be a 
method more effective than Method F that does not rely on the closed loop 
matrix being diagonalizable. Other possible research lines are applying these 
approaches to discrete-time Riccati equations (DARE) or more generally to 
non-symmetric algebraic Riccati equations (NARE). 
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